% replication_code_figure_all.m
% 
% This matlab code replicates figures 1-4 in Miyamoto, W., T.L. Nguyen, 
% and H. Oh (2023), "In Search of Dominant Drivers of the Real Exchange 
% Rate," Review of Economics and Statistics.
% 




%% Setup code to plot Figures 1-2
close all; clear all;
addpath('mat_files');
addpath('../Function'); 
% "mat_files": folder where mat files generated from this code are stored
% "../Function": folder with scripts/functions used in this code

setup_workspace;

s_LineWidth  = 3;
s_LineWidth2 = 2;
s_LineWidth3 = 2;
s_LineWidth4 = 1.5;
s_LineWidth5 = 1.5;

s_FontSize1  = 15;
s_FontSize2  = 12;
s_FontSize3  = 18;
s_FontSize4  = 20;

s_color1     = [0.0 0.0 0.0];
s_color2     = [1.0 0.0 0.0];
s_color3     = [0.3 0.4 0.7];
s_color4     = [0.3 0.3 0.3];
s_color5     = [0.7 0.0 0.0];

face1        = 0.9*[1 1 1];
face3        = 0.3*[1 1 1];

% percentile range
i_up   = 84; 
i_down = 16;


%% load files for VAR

jj_case=[814:820] ; % case to run VAR
countryname = {'USA','CAN','JPN','GBR','DEU','FRA','ITA'};

% start iteration
count=0;
for jj=jj_case
    
    count=count+1;    
    
    % load VAR
    ps_tmp.fpath.cucd=pwd; % current folder
    ps_tmp.file_name=['var_case' num2str(jj)];
    cd([ps.fpath.Case '/' ps_tmp.file_name]);
    load([ps_tmp.file_name '_var_result'])
    cd(ps_tmp.fpath.cucd)

    setup_workspace;
    
    % obtain impulse response
    opt.T_irf=20+1+3; % IRF length
    opt.scale_irf=1; % scale of IRF
    
    % case of MBC identification: k-th variable
    for k_th=1:ps.N % k-th variable

        Y_irf_draws=zeros(opt.T_irf,ps.N,ps.nsave); % preallocation
        % horizon by variable by draw

        for i=1:ps.nsave % draws
            % companion form coefficient hx : X(t)=hx*X(t-1)+netashock*e(t)
            hx=zeros(ps.N*ps.p,ps.N*ps.p); % preallocation
            hx(ps.N+1:end,1:ps.N*(ps.p-1))=eye((ps.p-1)*ps.N);

            tmp=out.ALPHA_draws(2:end,:,i); % take out constants
            hx(1:ps.N,1:ps.N*ps.p)=transpose(tmp); % transpose 

            %initial condition
            % x0=zeros(ps.N*ps.p,1);
            % % Cholesky decomposition
            % x0(1:ps.N,1)=chol(out.SIGMA_draws(:,:,i),'lower')*rs.q_draws(:,i);

            % create shock coefficient netashock 
            netashock=zeros(ps.N*ps.p,ps.N);
            netashock(1:ps.N,1)=chol(out.SIGMA_draws(:,:,i),'lower')*rs.q_draws{k_th}(:,i); % first shock

            % set control variables 
            gx=eye(ps.N,ps.N*ps.p);

            % create shock : shock at time t=2
            e=zeros(opt.T_irf,ps.N); % time by shock
            e(2,1)=1*opt.scale_irf; % 1st shock = MBC shock

            % compute irf
            [Y,~]=simu_1st(gx,hx,netashock,1,e);
            Y_irf_draws(:,:,i)=Y; % horizon by variable by draw
        end

        % keep irf
        Y_irf_draws_keep1{count,k_th}=Y_irf_draws(2:end,:,:); 
        Y_irf_draws(:,5,:)=-Y_irf_draws(:,5,:); % flip rer
        Y_irf_draws_keep{count,k_th}=Y_irf_draws(2:end,:,:); 
        % {case, k_th MBC shock}(horizon, variable, draw)
    end
end

%% VAR variable ordering (also holds for MBC shock ordering)
% 1. Y (output)
% 2. C (consumption)
% 3. h (hours worked)
% 4. NXY (net-export-to-output ratio)
% 5. RER (real exchange rate)
% 6. PI (inflation)
% 7. INT (nominal interest rate)
% 8. I (capital investment)

%--------------------------------------------------------------------------
%% Calculate the conditional UIP wedges
% UIP_irf = nom. int. rate diff. + \Delta RER(+1) - inflation diff.
% UIP_irf_draws(2:end,i) = Y(2:opt.T_irf-1,7) + ( Y(3:opt.T_irf,5) - Y(2:opt.T_irf-1,5)) - Y(3:opt.T_irf,6);

%--------------------------------------------------------------------------
% 1. UIP wedge conditional on dominant output shock
UIP_irf_draws_Y = NaN(opt.T_irf,ps.nsave,7);
for i=1:7
    data = squeeze(Y_irf_draws_keep1{i,1});
    UIP_irf_draws_Y(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
end
%--------------------------------------------------------------------------
% 2. UIP wedge conditional on dominant RER shock
UIP_irf_draws = zeros(opt.T_irf,ps.nsave,7);
for i=1:7
    data = squeeze(Y_irf_draws_keep1{i,5});
    UIP_irf_draws(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
end
%save('UIP_wedge_RERshock','UIP_irf_draws')% Save UIP wedge for RER shock for future reference
%--------------------------------------------------------------------------
% 3. UIP wedge conditional onl dominant INT shock
UIP_irf_draws_INT = zeros(opt.T_irf,ps.nsave,7);
for i=1:7
    data=squeeze(Y_irf_draws_keep{i,7});
    UIP_irf_draws_INT(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
end
%--------------------------------------------------------------------------
% 4. UIP wedge conditional on dominant PI shock
UIP_irf_draws_PI = zeros(opt.T_irf,ps.nsave,7);
for i=1:7
    data = squeeze(Y_irf_draws_keep1{i,6});
    UIP_irf_draws_PI(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
end
%--------------------------------------------------------------------------
% 5. UIP wedge conditional on dominant C shock
UIP_irf_draws_C = zeros(opt.T_irf,ps.nsave,7);
for i=1:7
    data = squeeze(Y_irf_draws_keep1{i,2});
    UIP_irf_draws_C(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
end
%--------------------------------------------------------------------------
% 6. UIP wedge conditional on dominant h shock
UIP_irf_draws_H = zeros(opt.T_irf,ps.nsave,7);
for i=1:7
    data = squeeze(Y_irf_draws_keep1{i,3});
    UIP_irf_draws_H(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
end
%--------------------------------------------------------------------------
% 7. UIP wedge conditional on dominant I shock
UIP_irf_draws_I = zeros(opt.T_irf,ps.nsave,7);
for i=1:7
    data = squeeze(Y_irf_draws_keep1{i,8});
    UIP_irf_draws_I(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
end
%--------------------------------------------------------------------------
% 8. UIP wedge conditional on dominant NXY shock
UIP_irf_draws_NXY = zeros(opt.T_irf,ps.nsave,7);
for i=1:7
    data = squeeze(Y_irf_draws_keep1{i,4});
    UIP_irf_draws_NXY(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
end


%% change the first impact 
opt.T_irf=20+1; %20+1; % IRF length


%% Plot Figure 1A
close all

figure1=figure();
s_col=5;
s_row=2;
s_time=0:opt.T_irf-1; % period for irf
         
i_var_all = [1 2 3 8 4]; % response variable index (Y,C,h,I,NXY)
i_MBC_all = [1 2 3 8];   % target variable index (Y,C,h,I)     
titlename={'Rel. Output','Rel. Consumption','Rel. Hours Worked','Rel. Investment','Net Exports','Real Exchange Rate'};

% USA
kk=1; 
for ii=1:5
    i_var=i_var_all(ii);
            
    subplot(s_row,s_col,ii)
    high = prctile(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),i_up,2);
    low  = prctile(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),i_down,2);
    h = fill([s_time fliplr(s_time)],[high' fliplr(low')],face1);
    set(h,'EdgeColor',face3)
    %set(h,'LineStyle','-.')
    hold on

    hp2 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),2),'-','Color',s_color1,'LineWidth',s_LineWidth);
    hp3 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(2)}(s_time+1,i_var,:)),2),'-.','Color',s_color4,'LineWidth',s_LineWidth);
    hp4 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(3)}(s_time+1,i_var,:)),2),'--','Color',s_color4,'LineWidth',s_LineWidth);
    hp5 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(4)}(s_time+1,i_var,:)),2),'*-','Color',s_color4,'LineWidth',s_LineWidth);
    hold off
    
    % set ylim
    if ii==1 || ii==2 || ii==3
        ylim([-0.2 1.0])
    elseif ii==4
        ylim([-0.5 2.2])
    elseif ii==5
        ylim([-0.2 0.05])
    end
            
    line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)
           
    xlabel('Quarter','FontSize',s_FontSize1)   
    if kk==1
        title(titlename{ii},'FontSize',s_FontSize3) 
    end
    if kk==1
        if ii==1
            legend([hp2 hp3 hp4 hp5], 'Y shock','C shock','h shock','I shock')   %'NXY shock','REER shock')
            set(legend,'FontSize',s_FontSize1,'Location','Northeast','color','none')
            legend boxoff
         end
    end
    if ii==1
        ylabel(countryname{kk},'FontSize',s_FontSize4,'fontweight','bold')     
    end
                        
    set(gca,'FontSize',s_FontSize3)
end 
     
% Median G6
for ii=1:5
    i_var=i_var_all(ii);
    
    for kk=2:7 
        mediang61(:,kk-1) = median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),2);
        mediang62(:,kk-1) = median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(2)}(s_time+1,i_var,:)),2);
        mediang63(:,kk-1) = median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(3)}(s_time+1,i_var,:)),2);
        mediang64(:,kk-1) = median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(4)}(s_time+1,i_var,:)),2);
    end
    
    subplot(s_row,s_col,ii+5)
    hold on
    hp2 = plot(s_time,median(mediang61,2),'-','Color',s_color1,'LineWidth',s_LineWidth);
    hp3 = plot(s_time,median(mediang62,2),'-.','Color',s_color4,'LineWidth',s_LineWidth);
    hp4 = plot(s_time,median(mediang63,2),'--','Color',s_color4,'LineWidth',s_LineWidth);
    hp5 = plot(s_time,median(mediang64,2),'*-','Color',s_color4,'LineWidth',s_LineWidth);
    
    % set ylim
    if ii==1 || ii==2 || ii==3
        ylim([-0.2 1.0])
    elseif ii==4
        ylim([-0.5 2.2])
    elseif ii==5
        ylim([-0.2 0.05])
    end
    
    line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)
    hold off

    xlabel('Quarter','FontSize',s_FontSize1)   
    if kk==1
        title(titlename{ii},'FontSize',s_FontSize3) 
    end
    if kk==1
        if ii==1
            legend([hp2 hp3 hp4 hp5], 'Y shock','C shock','h shock','I shock')   %'NXY shock','REER shock')
            set(legend,'FontSize',s_FontSize1,'Location','Northeast','color','none')
            legend boxoff
        end
    end
    if ii==1
        ylabel('G6 Median','FontSize',s_FontSize4,'fontweight','bold')     
    end
    set(gca,'Box','on');             
    set(gca,'FontSize',s_FontSize3)
end

set(gcf, 'PaperPosition', [0 0 18*5/4 8]); %Position plot at left hand corner with width 5 and height 5.
set(gcf, 'PaperSize', [18*5/4 8]);         %Set the paper to have width 5 and height 5.       

% store figures in "figures" folder
cd('figures')
name1=strcat(['fig_1a']);
print(name1, '-dpdf')
cd('..')


%% Plot Figure 1B: For USA and median G6, plot irf for PI, RER, INT, UIP (Y, C, h, I shocks)
close all

figure1=figure();
s_col=4;
s_row=2;
s_time=0:opt.T_irf-1; % period for irf

         
i_var_all = [6 5 7]; % response variable index (PI,RER,INT) + UIP (later)
i_MBC_all = [1 2 3 8];   % target variable index (Y,C,h,I)     
titlename={'Rel. Inflation Rate','Real Exchange Rate','Rel. Interest Rate','UIP Wedge'};
         
% USA
kk=1; 
for ii=1:3
    i_var=i_var_all(ii);
            
    subplot(s_row,s_col,ii)
    high = prctile(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),i_up,2);
    low  = prctile(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),i_down,2);
    h = fill([s_time fliplr(s_time)],[high' fliplr(low')],face1);
    set(h,'EdgeColor',face3)
    %set(h,'LineStyle','-.')
    hold on

    hp2 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),2),'-','Color',s_color1,'LineWidth',s_LineWidth);
    hp3 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(2)}(s_time+1,i_var,:)),2),'-.','Color',s_color4,'LineWidth',s_LineWidth);
    hp4 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(3)}(s_time+1,i_var,:)),2),'--','Color',s_color4,'LineWidth',s_LineWidth);
    hp5 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(4)}(s_time+1,i_var,:)),2),'*-','Color',s_color4,'LineWidth',s_LineWidth);
    hold off
    
    % set ylim
    if ii==1
        ylim([-0.2 0.1])
    elseif ii==2
        ylim([-1 1])
    elseif ii==3
        ylim([-0.02 0.1])
    end
            
    line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)
           
    xlabel('Quarter','FontSize',s_FontSize1)   
    if kk==1
        title(titlename{ii},'FontSize',s_FontSize3) 
    end
    if kk==1
        if ii==1
            legend([hp2 hp3 hp4 hp5], 'Y shock','C shock','h shock','I shock')   %'NXY shock','REER shock')
            set(legend,'FontSize',s_FontSize1,'Location','Southeast','color','none')
            legend boxoff
         end
    end
    if ii==1
        ylabel(countryname{kk},'FontSize',s_FontSize4,'fontweight','bold')     
    end
    set(gca,'FontSize',s_FontSize3)
end

ii=4; % UIP wedge
subplot(s_row,s_col,ii)
high = prctile(squeeze(UIP_irf_draws_Y(s_time+1,:,kk)),i_up,2);
low  = prctile(squeeze(UIP_irf_draws_Y(s_time+1,:,kk)),i_down,2);
h = fill([s_time fliplr(s_time)],[high' fliplr(low')],face1);
set(h,'EdgeColor',face3)

hold on
plot(s_time,squeeze(median(UIP_irf_draws_Y(s_time+1,:,kk),2)),'-','Color',s_color1,'LineWidth',s_LineWidth)
plot(s_time,squeeze(median(UIP_irf_draws_C(s_time+1,:,kk),2)),'-.','Color',s_color4,'LineWidth',s_LineWidth)
plot(s_time,squeeze(median(UIP_irf_draws_H(s_time+1,:,kk),2)),'--','Color',s_color4,'LineWidth',s_LineWidth)
plot(s_time,squeeze(median(UIP_irf_draws_I(s_time+1,:,kk),2)),'*-','Color',s_color4,'LineWidth',s_LineWidth)
hold off
    
line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)
   
ylim([-0.2 0.5])

xlabel('Quarter','FontSize',s_FontSize1)   
set(gca,'FontSize',s_FontSize3)
 
if kk==1
    title(titlename{ii},'FontSize',s_FontSize3) 
end

% Median G6
for ii=1:3
    i_var=i_var_all(ii);
    
    for kk=2:7 
        mediang61(:,kk-1) = median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),2);
        mediang62(:,kk-1) = median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(2)}(s_time+1,i_var,:)),2);
        mediang63(:,kk-1) = median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(3)}(s_time+1,i_var,:)),2);
        mediang64(:,kk-1) = median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(4)}(s_time+1,i_var,:)),2);
    end
    
    subplot(s_row,s_col,ii+4)
    hold on
    hp2 = plot(s_time,median(mediang61,2),'-','Color',s_color1,'LineWidth',s_LineWidth);
    hp3 = plot(s_time,median(mediang62,2),'-.','Color',s_color4,'LineWidth',s_LineWidth);
    hp4 = plot(s_time,median(mediang63,2),'--','Color',s_color4,'LineWidth',s_LineWidth);
    hp5 = plot(s_time,median(mediang64,2),'*-','Color',s_color4,'LineWidth',s_LineWidth);
    
    % set ylim
    if ii==1
        ylim([-0.2 0.1])
    elseif ii==2
        ylim([-1 1])
    elseif ii==3
        ylim([-0.02 0.1])
    end
    
    line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)
    hold off

    xlabel('Quarter','FontSize',s_FontSize1)   
    if kk==1
        title(titlename{ii},'FontSize',s_FontSize3) 
    end
    if kk==1
        if ii==1
            legend([hp2 hp3 hp4 hp5], 'Y shock','C shock','h shock','I shock')   %'NXY shock','REER shock')
            set(legend,'FontSize',s_FontSize1,'Location','Southeast','color','none')
            legend boxoff
        end
    end
    if ii==1
        ylabel('G6 Median','FontSize',s_FontSize4,'fontweight','bold')     
    end
    set(gca,'Box','on');
             
    set(gca,'FontSize',s_FontSize3)
end

ii=4; % UIP wedge

subplot(s_row,s_col,ii+4)
hold on
plot(s_time,median(squeeze(median(UIP_irf_draws_Y(s_time+1,:,2:7),2)),2),'-','Color',s_color1,'LineWidth',s_LineWidth)
plot(s_time,median(squeeze(median(UIP_irf_draws_C(s_time+1,:,2:7),2)),2),'-.','Color',s_color4,'LineWidth',s_LineWidth)
plot(s_time,median(squeeze(median(UIP_irf_draws_H(s_time+1,:,2:7),2)),2),'--','Color',s_color4,'LineWidth',s_LineWidth)
plot(s_time,median(squeeze(median(UIP_irf_draws_I(s_time+1,:,2:7),2)),2),'*-','Color',s_color4,'LineWidth',s_LineWidth)

line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)
hold off

ylim([-0.2 0.5])

xlabel('Quarter','FontSize',s_FontSize1)   

set(gca,'Box','on');
             
set(gca,'FontSize',s_FontSize3)

set(gcf, 'PaperPosition', [0 0 18*4/4 8]); %Position plot at left hand corner with width 5 and height 5.
set(gcf, 'PaperSize', [18*4/4 8]);         %Set the paper to have width 5 and height 5.       

% store figures in "figures" folder
cd('figures')
name1=strcat(['fig_1b']);
print(name1, '-dpdf')
cd('..')


%% Plot Figure 2
close all

figure1=figure();
s_col=4;
s_row=2;
s_time=0:opt.T_irf-1; % period for irf

i_var_all = [5 6 7 8 1 2 3 4]; % response variable index (RER, PI, INT, UIP wedge, Y, C, H, NXY)
i_MBC_all = [5];            % target variable index (RER)      
titlename = {'Real Exchange Rate','Rel. Inflation Rate','Rel. Interest Rate','UIP Wedge','Rel. Output','Rel. Consumption','Rel. Hours Worked','Net Exports'};
         
% USA
kk=1;
for ii=[1:3 5:8]
    i_var=i_var_all(ii);
    
    subplot(s_row,s_col,ii)
    high = prctile(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),i_up,2);
    low  = prctile(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),i_down,2);
    h = fill([s_time fliplr(s_time)],[high' fliplr(low')],face1);
    set(h,'EdgeColor',face3)
    %set(h,'LineStyle','-.')
    
    hold on
    hp1 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),2),'-','Color',s_color1,'LineWidth',s_LineWidth);
    % Other G6      
    for kkk=2:7
        mediang6(:,kkk-1)=median(squeeze(Y_irf_draws_keep{kkk,i_MBC_all(1)}(s_time+1,i_var,:)),2);
    end
    % Plot median G6
    hp2 = plot(s_time,median(mediang6,2),'-+','Color',s_color4,'LineWidth',s_LineWidth);
    hold off
        
    xlabel('Quarter','FontSize',s_FontSize1)   
    
    line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)
    
    if ii==1 || ii==5
        ylabel('Percent','FontSize',s_FontSize4)     
    end
    if ii==1
        legend([hp1 hp2 ], 'USA','G6 Median')   %'NXY shock','REER shock')
        set(legend,'FontSize',s_FontSize1,'Location','Southeast','color','none')
        legend boxoff
    end
             
    if i_var==1 || i_var==2 || i_var==3
        ylim([-0.2 0.6])
    end
    set(gca,'FontSize',s_FontSize3)
    if kk==1
        title(titlename{ii},'FontSize',s_FontSize3) 
    end
end

% UIP wedge
ii=4;

subplot(s_row,s_col,ii)
high = prctile(squeeze(UIP_irf_draws(s_time+1,:,kk)),i_up,2);
low  = prctile(squeeze(UIP_irf_draws(s_time+1,:,kk)),i_down,2);
h = fill([s_time fliplr(s_time)],[high' fliplr(low')],face1);
set(h,'EdgeColor',face3)

hold on
plot(s_time,squeeze(median(UIP_irf_draws(s_time+1,:,kk),2)),'-','Color',s_color1,'LineWidth',s_LineWidth)
plot(s_time,median(squeeze(median(UIP_irf_draws(s_time+1,:,2:7),2)),2),'-+','Color',s_color4,'LineWidth',s_LineWidth)
hold off

line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)

xlabel('Quarter','FontSize',s_FontSize1)   
set(gca,'FontSize',s_FontSize3)
 
if kk==1
    title(titlename{ii},'FontSize',s_FontSize3) 
end

set(gcf, 'PaperPosition', [0 0 18*4/4 8]); %Position plot at left hand corner with width 5 and height 5.
set(gcf, 'PaperSize', [18*4/4 8]);         %Set the paper to have width 5 and height 5.       

% store figures in "figures" folder
cd('figures')
name1=strcat(['fig_2']);
print(name1, '-dpdf')
cd('..')    



%% Setup code to plot Figure 3
setup_workspace

jj_case=[990:996]; % For Figure 3, we only use case 990 (USA)

count=0;
for jj=jj_case
    count=count+1;
    
    % load VAR estimation result files
    ps_tmp.fpath.cucd=pwd; % current folder
    ps_tmp.file_name=['var_case' num2str(jj)];
    cd([ps.fpath.Case '/' ps_tmp.file_name]);
    load([ps_tmp.file_name '_var_result'])
    cd(ps_tmp.fpath.cucd)

    setup_workspace;

    % impulse response
    opt.T_irf=20+1+3; %20+1; % IRF length
    opt.scale_irf=1; % scale of IRF
    
    % case of MBC identification : k-th variable
    for k_th=1:ps.N % k-th variable
        Y_irf_draws  = zeros(opt.T_irf,ps.N,ps.nsave); % preallocation
        Y_irf_draws2 = zeros(opt.T_irf,ps.N,ps.nsave); % preallocation
        Y_irf_draws3 = zeros(opt.T_irf,ps.N,ps.nsave); % preallocation
        % horizon by variable by draw

        for i=1:ps.nsave % draws
            % companion form coefficient hx : X(t)=hx*X(t-1)+netashock*e(t)
            hx=zeros(ps.N*ps.p,ps.N*ps.p); % preallocation
            hx(ps.N+1:end,1:ps.N*(ps.p-1))=eye((ps.p-1)*ps.N);

            tmp=out.ALPHA_draws(2:end,:,i); % take out constants
            hx(1:ps.N,1:ps.N*ps.p)=transpose(tmp); % transpose 

            %initial condition
            % x0=zeros(ps.N*ps.p,1);
            % % Cholesky decomposition
            % x0(1:ps.N,1)=chol(out.SIGMA_draws(:,:,i),'lower')*rs.q_draws(:,i);

            % create shock coefficient netashock 
            netashock=zeros(ps.N*ps.p,ps.N);
            netashock(1:ps.N,1)=chol(out.SIGMA_draws(:,:,i),'lower')*rs.q_draws{k_th}(:,i); % dominant output shock
            netashock2=zeros(ps.N*ps.p,ps.N);
            netashock2(1:ps.N,1)=chol(out.SIGMA_draws(:,:,i),'lower')*rs.q2_draws{k_th}(:,i); % orthogonalized dominant RER shock
            netashock3=zeros(ps.N*ps.p,ps.N);
            netashock3(1:ps.N,1)=chol(out.SIGMA_draws(:,:,i),'lower')*rs.q1_draws{k_th}(:,i); % dominant RER shock

            % set control variables 
            gx=eye(ps.N,ps.N*ps.p);

            % create shock : shock at time t=2
            e=zeros(opt.T_irf,ps.N); % time by shock
            e(2,1)=1*opt.scale_irf; % 1st shock = MBC shock

            % compute irf for three shocks
            [Y,~]=simu_1st(gx,hx,netashock,1,e);
            Y_irf_draws(:,:,i)=Y;
            [Y2,~]=simu_1st(gx,hx,netashock2,1,e);
            Y_irf_draws2(:,:,i)=Y2;
            [Y3,~]=simu_1st(gx,hx,netashock3,1,e);
            Y_irf_draws3(:,:,i)=Y3;
        end
        
        % dominant output shock
        Y_irf_draws_keep1{count,k_th}=Y_irf_draws(2:end,:,:); 
        Y_irf_draws(:,5,:)=-Y_irf_draws(:,5,:); % flip rer
        Y_irf_draws_keep{count,k_th}=Y_irf_draws(2:end,:,:); 
        
        % orthogonalized dominant RER shock
        Y_irf_draws_keep21{count,k_th}=Y_irf_draws2(2:end,:,:); 
        Y_irf_draws2(:,5,:)=-Y_irf_draws2(:,5,:); % flip rer
        Y_irf_draws_keep2{count,k_th}=Y_irf_draws2(2:end,:,:);
        
        % dominant RER shock
        Y_irf_draws_keep31{count,k_th}=Y_irf_draws3(2:end,:,:); 
        Y_irf_draws3(:,5,:)=-Y_irf_draws3(:,5,:); % flip rer
        Y_irf_draws_keep3{count,k_th}=Y_irf_draws3(2:end,:,:); 
    end
end

% Calculate uip wedge for G7 countries
UIP_irf_draws=zeros(opt.T_irf,ps.nsave,7);
UIP_irf_draws2=zeros(opt.T_irf,ps.nsave,7);
UIP_irf_draws3=zeros(opt.T_irf,ps.nsave,7);

for i=1:7
    data=squeeze(Y_irf_draws_keep1{i,5});
    UIP_irf_draws(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
    data=squeeze(Y_irf_draws_keep21{i,5});
    UIP_irf_draws2(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
    data=squeeze(Y_irf_draws_keep31{i,5});
    UIP_irf_draws3(1:end-2,:,i) = squeeze(data(1:opt.T_irf-2,7,:) + ( data(2:opt.T_irf-1,5,:) - data(1:opt.T_irf-2,5,:)) - data(2:opt.T_irf-1,6,:));
end

opt.T_irf=21;


%% Plot Figure 3
close all

figure1=figure();
s_col=4;
s_row=1;
s_time=0:opt.T_irf-1; % period for irf

i_var_all = [1 5 6 7 4]; % response variable index (Y, RER, PI, INT, NXY)
i_MBC_all = [5];         % target variable index (RER)
titlename = {'Rel. Output','Real Exchange Rate','Rel. Inflation Rate','Rel. Interest Rate','Net Exports','UIP Wedge'};

% USA       
kk=1;
for ii=1:4
    i_var=i_var_all(ii);
    
    subplot(s_row,s_col,ii)
    high = prctile(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),i_up,2);
    low  = prctile(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),i_down,2);
    h = fill([s_time fliplr(s_time)],[high' fliplr(low')],face1);
    set(h,'EdgeColor',face3)
    %set(h,'LineStyle','-.')
    
    hold on
    hp1 = plot(s_time,median(squeeze(Y_irf_draws_keep{kk,i_MBC_all(1)}(s_time+1,i_var,:)),2),'-.','Color',s_color1,'LineWidth',s_LineWidth);
    hp2 = plot(s_time,median(squeeze(Y_irf_draws_keep2{kk,i_MBC_all(1)}(s_time+1,i_var,:)),2),'-o','Color',s_color4,'LineWidth',s_LineWidth);
    hp3 = plot(s_time,median(squeeze(Y_irf_draws_keep3{kk,i_MBC_all(1)}(s_time+1,i_var,:)),2)+0.001,'-','Color',s_color4,'LineWidth',s_LineWidth);
    hold off
        
    xlabel('Quarter','FontSize',s_FontSize1)   
    
    line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)
    
    if ii==1
        ylabel('Percent','FontSize',s_FontSize4)     
    end
    
    set(gca,'FontSize',s_FontSize3)
    if kk==1
        title(titlename{ii},'FontSize',s_FontSize3) 
    end
    if ii==1
        ylim([-0.2 1.2])
        legend([hp1 hp2 hp3], 'Rel. Output Shock','Orthog. RER Shock','RER Shock')   %'NXY shock','REER shock')
        set(legend,'FontSize',s_FontSize1,'Location','Northeast','color','none')
        legend boxoff
    end
end

set(gcf, 'PaperPosition', [0 0 19.5*4/4 4]); %Position plot at left hand corner with width 5 and height 5.
set(gcf, 'PaperSize', [18.5*4/4 4]);         %Set the paper to have width 5 and height 5.       

% store figures in "figures" folder
cd('figures')
name1=strcat(['fig_3']);
print(name1, '-dpdf')
cd('..')



%% Setup code to plot Figure 4
%--------------------------------------------------------------------------
clear all;close all;clc
setup_workspace;
setup_seed;


%% Load model impulse responses
% load data from dynare
cd(ps.fpath.Main);
cd(ps.fpath.Models)
load('gnk2_loop4_irf_results.mat');
cd(ps.fpath.Main);

%%% plot dynare
close all; clear irf

% choose variables
var_name={'a1','a2','y1','y2','c1',...
       'c2','h1','h2','I1','I2',...
      'nxy','qD','pi1','pi2','r1',...
      'r2','wedge','psi1','mo1','mo2'};

% choose shock
%%%var_shock=['e_a1']; % productivity
%%%var_shock=['e_mo1']; % monetary
var_shock=['e_psi1']; % financial

for i=1:length(var_name)
    eval(['irf(i,:)=oo_.irfs.' char(var_name{i}) '_' char(var_shock) ';'])
end

% add 0 initial
irf=[zeros(size(irf,1),1) irf];
s_period=1:21+1;

% get model IRFs
IRFmodel_Y     = - (irf(3,s_period)-irf(4,s_period)); % store Y
IRFmodel_qD    = - irf(12,s_period); % store qD
IRFmodel_PI    = - (irf(13,s_period)-irf(14,s_period)); % store inflation
IRFmodel_INT   = -(irf(15,s_period)-irf(16,s_period)); % store interest rate
IRFmodel_wedge = - irf(17,s_period); % store UIP wedge
IRFmodel_NXY   = - irf(11,s_period);

%% Load model simulated data
load('all_workspace.mat')

setup_workspace;

s_case=6; % case of the dynare model
s_mbc=5;  % target variable for mbc

clear tmp
for i=1:ps.nsim
    tmp(:,:,i)=Y_irf_draws_all{s_case,i,s_mbc}; % case by simulation by mbc variable    
end

% percentile
tmp50 = prctile(tmp,50,3); tmp50 = tmp50';
tmp10 = prctile(tmp,10,3); tmp10 = tmp10';
tmp90 = prctile(tmp,90,3); tmp90 = tmp90';
tmp16 = prctile(tmp,16,3); tmp16 = tmp16';
tmp84 = prctile(tmp,84,3); tmp84 = tmp84';

% Y
IRFsim_Y_50 = tmp50(1,:); IRFsim_Y_16 = tmp16(1,:); IRFsim_Y_84 = tmp84(1,:);
% C
IRFsim_C_50 = tmp50(2,:); IRFsim_C_16 = tmp16(2,:); IRFsim_C_84 = tmp84(2,:);
% H
IRFsim_H_50 = tmp50(3,:); IRFsim_H_16 = tmp16(3,:); IRFsim_H_84 = tmp84(3,:);
% NXY
IRFsim_NXY_50 = tmp50(4,:); IRFsim_NXY_16 = tmp16(4,:); IRFsim_NXY_84 = tmp84(4,:);
% qD
IRFsim_qD_50 = -tmp50(5,:); IRFsim_qD_16 = -tmp16(5,:); IRFsim_qD_84 = -tmp84(5,:);
% Inflation
IRFsim_PI_50 = tmp50(6,:); IRFsim_PI_16 = tmp16(6,:); IRFsim_PI_84 = tmp84(6,:);
% Interest rate
IRFsim_INT_50 = tmp50(7,:); IRFsim_INT_16 = tmp16(7,:); IRFsim_INT_84 = tmp84(7,:);
% UIP wedge
IRFsim_wedge_50 = tmp50(9,:); IRFsim_wedge_16 = tmp16(9,:); IRFsim_wedge_84 = tmp84(9,:);
% Investment
IRFsim_I_50 = tmp50(8,:); IRFsim_I_16 = tmp16(8,:); IRFsim_I_84 = tmp84(8,:);


%% load empirical data
load('fig_irf_medianG7_RER_data.mat')

%% scale to match the impact RER response to a dominant RER shock
scale = IRFmodel_qD(2)/IRFsim_qD_50(2); % scale to match impact RER response
scale = scale/4; scale_model = 1/4;

%% Plot Figure 4
close all

s_LineWidth=3;
s_LineWidth2=2;
s_LineWidth3=2;
s_LineWidth4=1.5;
s_LineWidth5=1.5;

s_FontSize1=15;
s_FontSize2=12;
s_FontSize3=18;
s_FontSize4=20;

s_color1=[0 0 0];
s_color2=[1 0 0];
s_color3=[0.30 0.40 0.70];
s_color4=0.3*[1 1 1];

face3=0.3*[1 1 1];
face1=0.9*[1 1 1];


figure1=figure();
close all
s_col=5;
s_row=1;
s_time=0:20; % period for irf
T_irf = length(s_time);

ii = 0;

% UIP wedge plot
ii=ii+1;
subplot(s_row,s_col,ii)
high = scale*IRFsim_wedge_84(2:T_irf+1);
low  = scale*IRFsim_wedge_16(2:T_irf+1);  
h=fill([s_time fliplr(s_time)],[high fliplr(low)],face1);
set(h,'EdgeColor',face3)
%set(h,'LineStyle','-.')

hold on
hp1 = plot(s_time,scale*IRFsim_wedge_50(2:T_irf+1),'-','Color',s_color1,'LineWidth',s_LineWidth);
% hp2 = plot(s_time,median(squeeze(Y_irf_draws_keep{i_case_all1{kk}(:,ii),i_MBC}(:,i_var,:)),2),'-.','Color',[0 0 1],'LineWidth',s_LineWidth);
hp3 = plot(s_time,scale_model*IRFmodel_wedge(2:T_irf+1),'--','Color',s_color4,'LineWidth',s_LineWidth);
hp4 = plot(s_time,median_g7_data(1:T_irf,ii),'-.','Color',s_color4,'LineWidth',s_LineWidth);

 ylim([-1 0.2])
xlabel('Quarter','FontSize',s_FontSize1)
ylabel('Percent','FontSize',s_FontSize1)

line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)

legend([hp1 hp3 hp4], 'Simulated Data','Financial Shock','G7 Median') 
set(legend,'FontSize',s_FontSize1,'Location','Southeast','color','none')
legend boxoff

set(gca,'FontSize',s_FontSize3)
title('UIP wedge','FontSize',s_FontSize3) 


% RER plot
ii=ii+1;
subplot(s_row,s_col,ii)
high = scale*IRFsim_qD_84(2:T_irf+1);
low  = scale*IRFsim_qD_16(2:T_irf+1);  
h=fill([s_time fliplr(s_time)],[high fliplr(low)],face1);
set(h,'EdgeColor',face3)
%set(h,'LineStyle','-.')

hold on
hp1 = plot(s_time,scale*IRFsim_qD_50(2:T_irf+1),'-','Color',s_color1,'LineWidth',s_LineWidth);
% hp2 = plot(s_time,median(squeeze(Y_irf_draws_keep{i_case_all1{kk}(:,ii),i_MBC}(:,i_var,:)),2),'-.','Color',[0 0 1],'LineWidth',s_LineWidth);
hp3 = plot(s_time,scale_model*IRFmodel_qD(2:T_irf+1),'--','Color',s_color4,'LineWidth',s_LineWidth);
hp4 = plot(s_time,median_g7_data(1:T_irf,ii),'-.','Color',s_color4,'LineWidth',s_LineWidth);

%ylim([0 1])
xlabel('Quarter','FontSize',s_FontSize1)
%ylabel('Percent','FontSize',s_FontSize1)

line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)


set(gca,'FontSize',s_FontSize3)
title('Real Exchange Rate','FontSize',s_FontSize3) 


% Inflation plot
ii=ii+1;
subplot(s_row,s_col,ii)
high = scale*IRFsim_PI_84(2:T_irf+1);
low  = scale*IRFsim_PI_16(2:T_irf+1);  
h=fill([s_time fliplr(s_time)],[high fliplr(low)],face1);
set(h,'EdgeColor',face3)
%set(h,'LineStyle','-.')

hold on
hp1 = plot(s_time,scale*IRFsim_PI_50(2:T_irf+1),'-','Color',s_color1,'LineWidth',s_LineWidth);
% hp2 = plot(s_time,median(squeeze(Y_irf_draws_keep{i_case_all1{kk}(:,ii),i_MBC}(:,i_var,:)),2),'-.','Color',[0 0 1],'LineWidth',s_LineWidth);
hp3 = plot(s_time,scale_model*IRFmodel_PI(2:T_irf+1),'--','Color',s_color4,'LineWidth',s_LineWidth);
hp4 = plot(s_time,median_g7_data(1:T_irf,ii),'-.','Color',s_color4,'LineWidth',s_LineWidth);

%ylim([0 1])
xlabel('Quarter','FontSize',s_FontSize1)
%ylabel('Percent','FontSize',s_FontSize1)

line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)


set(gca,'FontSize',s_FontSize3)
title('Rel. Inflation Rate','FontSize',s_FontSize3) 


% Interest rate plot
ii=ii+1;
subplot(s_row,s_col,ii)
high = scale*IRFsim_INT_84(2:T_irf+1);
low  = scale*IRFsim_INT_16(2:T_irf+1);  
h=fill([s_time fliplr(s_time)],[high fliplr(low)],face1);
set(h,'EdgeColor',face3)
%set(h,'LineStyle','-.')

hold on
hp1 = plot(s_time,scale*IRFsim_INT_50(2:T_irf+1),'-','Color',s_color1,'LineWidth',s_LineWidth);
% hp2 = plot(s_time,median(squeeze(Y_irf_draws_keep{i_case_all1{kk}(:,ii),i_MBC}(:,i_var,:)),2),'-.','Color',[0 0 1],'LineWidth',s_LineWidth);
hp3 = plot(s_time,scale_model*IRFmodel_INT(2:T_irf+1),'--','Color',s_color4,'LineWidth',s_LineWidth);
hp4 = plot(s_time,median_g7_data(1:T_irf,ii),'-.','Color',s_color4,'LineWidth',s_LineWidth);

%ylim([0 1])
xlabel('Quarter','FontSize',s_FontSize1)
%ylabel('Percent','FontSize',s_FontSize1)

line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)


set(gca,'FontSize',s_FontSize3)
title('Rel. Interest Rate','FontSize',s_FontSize3) 


% Net exports plot
ii=ii+1;
subplot(s_row,s_col,ii)
high = scale*IRFsim_NXY_84(2:T_irf+1);
low  = scale*IRFsim_NXY_16(2:T_irf+1);  
h=fill([s_time fliplr(s_time)],[high fliplr(low)],face1);
set(h,'EdgeColor',face3)
%set(h,'LineStyle','-.')

hold on
hp1 = plot(s_time,scale*IRFsim_NXY_50(2:T_irf+1),'-','Color',s_color1,'LineWidth',s_LineWidth);
% hp2 = plot(s_time,median(squeeze(Y_irf_draws_keep{i_case_all1{kk}(:,ii),i_MBC}(:,i_var,:)),2),'-.','Color',[0 0 1],'LineWidth',s_LineWidth);
hp3 = plot(s_time,scale_model*IRFmodel_NXY(2:T_irf+1),'--','Color',s_color4,'LineWidth',s_LineWidth);
hp4 = plot(s_time,median_g7_data(1:T_irf,ii),'-.','Color',s_color4,'LineWidth',s_LineWidth);

%ylim([0 1])
xlabel('Quarter','FontSize',s_FontSize1)
%ylabel('Percent','FontSize',s_FontSize1)

line([s_time(1) s_time(end)],[0 0],'Color',s_color4,'LineStyle',':','LineWidth',s_LineWidth5)

set(gca,'FontSize',s_FontSize3)
title('Net Exports','FontSize',s_FontSize3) 

set(gcf, 'PaperPosition', [0 0 18*5/4 4]); %Position plot at left hand corner with width 5 and height 5.
set(gcf, 'PaperSize', [18*5/4 4]);         %Set the paper to have width 5 and height 5.  


% store figures in "figures" folder
cd('figures')
name1=strcat(['fig_4']);
print(name1, '-dpdf')
cd('..')


